########################################################################################
########################################################################################
########  Does discrimination breed grievances--and do grievances breed violence?  #####
########            R Code for Main Figures                                        #####
########################################################################################
########################################################################################
rm(list=ls(all=TRUE))

library(foreign)
library(lme4)
library(scales)
library(ggplot2)
library(memisc)
library(Zelig)
library(MASS)


# SET PATHS AND LOAD DATA
setwd("~/Dropbox/GIGA/Religion_2 CMPS/")
data <- read.dta("ras_rcdc_long_130107_controls_grievances.dta")


# Figure 1: Historgram of the Discrimination Index 
qplot(mmx,data=data,geom=c("histogram"),ylab="Density",xlab="Discrimination Index",main="")+
  theme_bw()


#Substantive Effect
baseline <- glm(v95group ~ t1+ t2+ t3+ mz00pop+ reltype1+reltype2+reltype3+reltype4+reltype5+reltype6+
                  reltype7+reltype8+reltype9+reltype10+reltype11+reltype12+
                  reltype13+reltype14+reltype15+reltype16+reltype17+reltype18++zpoplog+ zpcgdpl+ zpolity2+ zdurable+ mmx_l1, family=binomial(link="logit"), data=data)
summary(baseline)

means <- coef(baseline)
var <- vcov(baseline)
draws <- mvrnorm(n=1000, means, var)
concvec <- seq(min(baseline$data$mmx_l1,na.rm=T), max(baseline$data$mmx_l1,na.rm=T), by=((max(baseline$data$mmx_l1,na.rm=T)-min(baseline$data$mmx_l1,na.rm=T))/(500-1)))

meanvec <- lowervec <- highervec <- NULL
for(i in 1:length(concvec)){
  out <- draws[,1] + 7*draws[,2]+49*draws[,3]+(7^3)*draws[,4]+mean(data$mz00pop,na.rm=T)*draws[,5]+
    1*draws[,18]+mean(data$zpoplog,na.rm=T)*draws[,24]+mean(data$zpcgdpl,na.rm=T)*draws[,25]+
    mean(data$zpolity2,na.rm=T)*draws[,26]+ mean(data$zdurable,na.rm=T)*draws[,27]+concvec[i]*draws[,28] 
  prob <- 1/(1+exp(-out))
  meanhelp <- means[1] + 7*means[2]+49*means[3]+(7^3)*means[4]+mean(data$mz00pop,na.rm=T)*means[5]+
    1*means[18]+mean(data$zpoplog,na.rm=T)*means[24]+mean(data$zpcgdpl,na.rm=T)*means[25]+
    mean(data$zpolity2,na.rm=T)*means[26]+ mean(data$zdurable,na.rm=T)*means[27]+concvec[i]*means[28]
  mean <- 1/(1+exp(-meanhelp))
  lower <- quantile(prob,0.05)
  higher <- quantile(prob,0.95)
  meanvec <- c(meanvec,mean)
  lowervec <- c(lowervec,lower)
  highervec <- c(highervec,higher)
}

pdf("DiscriminationEffect.pdf",height=7, width=7*1.62)
plot(concvec,meanvec,type="n",ylim=c(0,0.2), xlab="Discrimination", ylab="Probability of Grievance by Group")
polygon(c(concvec, rev(concvec)), c(lowervec,rev(highervec)), col="grey80", border=NA)
points(concvec,meanvec,type="l",lty=1)
rug(data$mmx_l1)
dev.off()









